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Abstract 

The  standard  harmonic  balance  method  consists  in  expanding  the  displacement  of  an  oscillator 
as  a  Fourier  cosine  series  in  time.  A  key  modification  is  proposed  here,  in  which  the  conservative 
force  is  additionally  expanded  as  a  Fourier  sine  series  in  space.  As  a  result,  the  steady-state 
oscillation  frequency  can  be  expressed  in  terms  of  a  Bessel  series,  and  the  sums  of  many  such 
series  are  known  or  can  be  developed.  The  method  is  illustrated  for  five  different  physical 
situations,  including  a  ball  rolling  inside  a  V-shaped  ramp,  an  electron  attracted  to  a  charged 
filament,  a  large-amplitude  pendulum,  and  a  Duffing  oscillator.  As  an  example  of  the  results, 
the  predicted  period  of  a  simple  pendulum  swinging  between  —90°  and  +90°  is  found  to  be 
only  0.4%  larger  than  the  exact  value.  Even  better,  the  predicted  frequency  for  the  V-ramp 
case  turns  out  to  be  exact. 

Keywords:  Fourier  expansion,  harmonic  balance,  Bessel  series,  work-energy  theorem,  nonlinear 
oscillator,  pendulum. 


1  Introduction 

Nonlinear  oscillators  are  ubiquitous  in  physical  and  engineering  systems.  Approximate  analytical 
solutions  are  useful  for  exploring  the  features  of  such  systems.  Simple  perturbation  theory,  however, 
can  give  rise  to  secular  terms  (i.e.,  ones  that  diverge  in  time)  which  result  in  unphysical  solutions. 
As  a  consequence,  alternative  approximations  are  needed.  One  of  these  is  the  Lindstedt-Poincare 
method  [I],  in  which  both  the  frequency  and  the  nonlinearity  are  expanded  as  power  series  in  an 
appropriate  parameter.  The  requirement  that  there  should  not  be  unbounded  terms  then  yields  a 
condition  for  the  amplitude-frequency  relationship.  J.-H.  He  has  developed  other  approaches  for 
analyzing  nonlinear  oscillators,  including  the  homotopy  and  the  variational  iteration  methods  [2]. 
More  widespread  is  the  harmonic  balance  method  [3],  wherein  the  displacement  is  expanded  as  a 
series  of  frequency  overtones.  The  equation  of  motion  imposes  a  relationship  between  the  amplitude 
and  the  frequency  of  the  system.  The  lowest  nonzero  term — the  harmonic — is  dominant. 
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In  the  present  article,  a  double  Fourier  method  is  developed  that  improves  upon  the  conventional 
harmonic  balance  approximation  [4,  5,  6].  Five  key  examples  of  physical  importance  are  presented 
here  to  illustrate  the  method  and  the  results  are  compared  with  solutions  obtained  using  a 
work-energy  principle.  In  the  course  of  applying  the  method  to  these  examples,  novel  sums  of 
certain  Bessel  series  are  developed  from  special  cases  tabulated  in  the  literature.  In  those  cases 
when  the  Bessel  series  can  be  summed  analytically,  the  new  method  in  effect  extends  standard 
harmonic  balance  to  include  all  overtones.  Even  when  the  Bessel  series  cannot  be  analytically 
summed,  the  first  few  terms  provide  insight  about  the  behavior  of  the  system. 


2  Approximate  frequency  using  double  Fourier  harmonic 
balance 

The  goal  of  this  section  is  to  estimate  the  frequency  of  oscillation,  an  important  parameter  of  the 
motion  [5] .  Write  the  equation  of  motion  for  the  displacement  r  as  a  function  of  time  t  in  the  form 


^  +  /(*)  =  ° 


(1) 


where  f{x)  is  the  negative  of  the  force  per  unit  mass  and  is  in  general  a  nonlinear  function  of  x. 
Assume  the  amplitude  of  oscillation  is  A,  with  the  displacement  varying  continuously  over  the  range 
— A  <  x  <  A.  Choosing  the  origin  so  that  /( 0)  =  0,  expand  f(x)  in  the  Fourier  sine  series 


/O)  =  ^2br 


sin 


/ nnx\ 

it) 


where  the  coefficients  are 


b™=‘jfo  /(sO  sin  (r^C)  dx. 

The  equation  of  motion  then  becomes 


a  x  . 


sin 


n—l 


(TH- 


(2) 


(3) 


(4) 


Assume  the  oscillator  starts  out  at  t  =  0  with  zero  initial  velocity  at  a :  =  A,  has  zero  displacement 
at  a  quarter  period  t  =  T/ 4,  just  attains  x  =  — A  with  zero  velocity  at  half  a  period,  crosses  back 
through  x  =  0  at  three-quarters  of  a  period,  and  the  cycle  starts  over  at  t  =  T.  In  that  case, 
the  displacement  can  be  written  as  a  cosine  series  in  time  with  only  odd  harmonics  due  to  the 
symmetry, 

x(t)  =  Aicoswf  +  A3  cos  Suit  +  A5  cos  5 uit  +  ...  (5) 

where  the  angular  frequency  of  oscillation  is  ui  =  2i r/T  and  the  total  amplitude  is 


A  —  A\  A3  +  A3  +  — 

Substitution  of  Eq.  (5)  into  (4)  leads  to 

E.  ,2  ■(A,  .  I"  (A  —  A3  —  ...)mr  cosiot  +  A3mrcos3ujt  +  ... 

Am(muj)  cosmojt  =  )  y  bn sin  - — - 

odd  m=l  n—l  L 


(6) 


(7) 
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where  Eq.  (6)  was  used  to  eliminate  A1  on  the  right-hand  side.  Defining  the  difference  between  the 
exact  and  linear-order  solution  as 


6(t)  =  E  (cos  mwt  —  cos  cot)  =  pp  —  cos  uit, 


odd  m— 3 


(8) 


Eq.  (7)  can  be  rewritten  as 

OO  OO 

ur  E  A2k+i(2k  +  l)2  [cos(2/e  +  1  )ut]  =  E  bnsin(mr  cos  uit  +  mrd).  (9) 

k= 0  n= 1 

This  expression  is  exact,  but  now  make  the  approximation  that  6  can  be  neglected  in  the  context 
of  Eq.  (9).  Expression  9.1.45  in  Abramowitz  and  Stegun  [7]  is 


sin  (zcosO)  =  2^  (— l)fc  J2k+i  ( z )  cos  [(2k  +  1)0]  (10) 

k= o 


in  which  J  is  a  Bessel  function  of  integer  order.  (This  expression  is  a  real-valued  Jacobi-Anger 
expansion.)  Consequently 


OO 

sin  (wr  cos  cot)  =  2  E  (-l)fc  J2k+i  (nn)  cos  [(2k  +  1  )w<] .  (11) 

k- o 

Substituting  this  series  into  the  right-hand  side  of  Eq.  (9)  with  <5  =  0  and  equating  coefficients  of 
the  time  harmonics  leads  to 


w2A2fc+ i  — 


2(— l)fe 
(2  k  +  l)2 


OO 

E  bnJ2k+i(mr). 

n= 1 


(12) 


Summing  both  sides  over  integer  values  of  k  from  0  to  oo  gives 


UJ 


2 


n= 1 


Ji(mr) 


1 

9 


J3(mr) 


1 

25 


J5(n7r)  -  ... 


(13) 


using  Eq.  (6).  Recalling  that  bn  is  determined  for  a  given  force  law  by  Eq.  (3),  this  expression 
predicts  the  angular  frequency  of  oscillation,  but  it  is  approximate  because  of  the  neglect  of  S  in 
Eq.  (9).  To  confirm  the  accuracy  of  this  approximation,  the  prediction  can  be  compared  to  the 
exact  frequency  for  various  physically  interesting  systems. 


3  Exact  frequency  calculated  using  the  work-kinetic-energy 
theorem 

The  acceleration  a  =  cPx/dt 2  is  related  to  the  velocity  v  =  dx/dt  by 


dv  dx  dv  dv 

dt  dt  dx  dx 
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and  therefore  Eq.  (1)  becomes 


vdv  =  —f(x)dx. 


Integrating  this  expression  from  v  =  0  at  x  =  A  to  velocity  v  at  displacement  x  gives 


=  /  f(x')dx'. 


(15) 


(16) 


The  left-hand  side  is  the  change  in  kinetic  energy  and  the  right-hand  side  is  the  work,  both  per 
unit  mass.  (For  clarity,  primes  have  been  added  to  the  dummy  variable  of  integration.)  Rewrite 
this  equation  as 


dx 

dt 


f(x')dx'. 


(17) 


Since  it  requires  one-quarter  of  a  period  T  for  the  oscillator  to  move  from  x  =  0  to  x  =  A, 


T 

~4 


dx 

]/2fxAf(x')dx' 


and  therefore  the  exact  angular  frequency  is 


k^exact  —  2  7T /T 


(18) 


(19) 


4  Physical  examples  applying  the  method 

4.1  The  simple  harmonic  oscillator 

Substituting  f(x)  =  ujgX  into  Eq.  (19)  and  changing  variables  to  9  according  to  x  =  A  cos  9,  the 
exact  angular  frequency  is  found  to  be  weXact  =  wo  as  expected.  Here,  for  example,  loq  =  k/m  for 
a  particle  of  mass  m  oscillating  on  a  spring  of  stiffness  k. 

Equation  (13)  becomes 


cc2  = 


k—0  K  J  n=  1 


X  Sill 


/ n7TX\ 

\~A~  ) 


dx 


using  Eq.  (3),  where  the  Fourier  series  is  that  of  a  sawtooth  wave,  and  thus 

( -i)k  ^  (-1  y" 


UJ 


=  4 w02  E 


E 


(2k  +  l)2 

fc= 0  v  ’  n—l 


~J2k+i(mr)- 


(20) 


(21) 


However,  Ec{.  (58)  in  Appendix  A  shows  that  only  the  k  =  0  term  is  nonzero,  leading  to 

OO  /_  1  \n+l 

w2  =  4wq  V - Ji(mr).  (22) 

mr 

n—l 
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Equation  (55)  in  Appendix  A  implies  that  this  sum  equals  |  and  consequently 


uj  —  ojq.  (23) 

A  simple  harmonic  oscillator  has  no  overtones  and  thus  A3  =  A5  =  ...  =  0,  so  that  5  =  0  according 
to  Eq.  (8).  As  a  result,  the  frequency  obtained  using  the  first-order  double  Fourier  harmonic 
balance  method  is  exact. 


4.2  Motion  down  two  intersecting  planes 


Consider  a  ball  starting  on  a  ramp  inclined  at  angle  <j>  relative  to  the  horizontal,  rolling  down  to  the 
origin,  and  then  rolling  up  another  ramp  inclined  at  the  same  angle  to  the  horizontal,  so  that  the 
two  ramps  form  a  V-shape.  In  the  absence  of  rolling  friction  or  air  drag,  the  ball  rolls  a  distance 
A  along  each  ramp.  Defining  ujq  =  gsmc/j/A,  where  g  is  the  gravitational  field,  the  restoring  force 
divided  by  the  mass  of  the  ball  is  the  step  function  f(x)  =  A£jgSgn(:r).  Inserting  this  force  law 
into  Eq.  (19)  and  changing  variables  to  6  according  to  x  =  A  cos2  0,  the  exact  angular  frequency  is 
found  to  be 


k^exact 


2^° 


l.lllWQ. 


(24) 


Since  ujq  depends  on  A,  as  will  also  be  the  case  for  every  remaining  example  treated  in  this  article, 
the  frequency  of  all  but  the  linear  oscillator  of  Sec.  4.1  in  general  depends  on  the  amplitude  of 
motion. 

Retaining  only  the  Ji  term,  Eqs.  (3)  and  (13)  predict 


U! 


A  9  00 


.  /  Tl'KX  \ 

sm  y  — —  j  ax 


(25) 


and  the  Fourier  integral  for  this  square  wave  gives 


UJ 


2 


odd  n— 1 


Ji(n7r) 

n 


(26) 


According  to  Eq.  (69)  in  Appendix  B,  the  sum  is  \  and  so  the  approximate  frequency  is 


1.128tdo 


(27) 


which  is  only  1.5%  larger  than  the  exact  answer.  By  retaining  all  of  the  terms  in  the  Bessel  series 
in  Eq.  (13)  and  using  Eq.  (74)  in  Appendix  C,  one  obtains 


UJ 


2 


4w| 

7T 


i  l 
25  '  5 


(28) 


Since  the  sum  of  the  alternating  reciprocals  of  the  odd  cubes  is  /3(3)  =  7t3/32,  where  /3  is  the 
Dirichlet  beta  function,  the  square  root  of  Eq.  (28)  is  equal  to  the  result  in  Eq.  (24).  This  result 
illustrates  that  in  the  present  example  S(t)  can  be  exactly  dropped  in  Eq.  (9)  even  though  it  is  not 
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Figure  1.  Graphs  of  the  net  displacement  from  Eq.  (29),  the  unit-amplitude 
lowest-order  Fourier  component  from  Eq.  (5),  and  the  difference  between  them  from 
Eq.  (8). 


zero  at  all  times.  Specifically,  S(t)  is  given  by  Eq.  (8)  as  graphed  in  Fig.  1,  where  the  displacement 
of  the  ball  is  given  by  the  slope-matched  parabolae 

x(t)  \i~^y  for2°^^f 

-f  =  {  -l  +  (f-2)2  for  ?<*<¥  (29) 

[  l-(f  -4)2  for^<t<T 

since  the  acceleration  down  the  ramps  has  the  constant  values  ±g  sin  <f>.  Thus  the  coefficients  of  the 
temporal  Fourier  series  in  Eq.  (5)  are 


1  -  32(~ir  fo0i 

A2m+1  “  7r3(2m  +  l)3  ’  (30) 

thereby  confirming  the  series  found  in  Eq.  (28)  above. 

4.3  Electric  force  inversely  proportional  to  distance 

An  electron  of  charge  —  e  and  mass  m  is  released  from  rest  a  distance  A  away  from  a  long  straight 
plasma  filament  of  positive  linear  charge  density  A.  Defining  uiq  =  2keX/mA2  where  k  is  the 
Coulomb  constant,  the  attractive  force  per  unit  mass  on  the  electron  is  f(x)  =  A2uiq/x.  Inserting 
this  force  law  into  Eq.  (19)  and  changing  variables  to  it  =  y/ln(A/x)  to  get  a  Gaussian  integral, 
the  exact  angular  frequency  is  found  to  be 


k^exact 


1.253u;o. 


(31) 
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To  apply  the  double  Fourier  harmonic  balance  method,  Eq.  (3)  becomes 


The  sine  integral  Si (t)  is  defined  as  the  integral  of  the  sine  function, 


so  that 

bn  =  2AugSi(n7r). 

Substituting  this  result  into  Eq.  (13),  the  approximate  frequency  is 


U! 


2 


k- 0 


(~l)fc 

(2k  +  l)2 


OO 

Y.  Si(mr)J2k+i(nn). 

n= 1 


(32) 


(33) 

(34) 


(35) 


The  partial  sums  over  n  oscillate  and  approach  !2  as  the  upper  limit  goes  to  infinity  for  any 
non-negative  integer  value  of  k.  For  example,  Fig.  2  plots  the  partial  sums  for  the  lowest  order  term 
corresponding  to  k  =  0, 

N 

S(N)  =  Y,Si(nn)Mnn),  (36) 

71=1 


Figure  2.  Graphs  of  the  partial  sums  in  Eq.  (36),  corresponding  to  k  =  0  in  Eq.  (35), 
for  values  of  N  ranging  from  one  to  a  million.  The  sums  approach  ^  from  above  for 
odd  values  of  N  and  from  below  for  even  values  of  N.  Similar  behavior  is  observed  for 
positive  even  values  of  k.  For  positive  odd  values  of  k,  the  partial  sums  for  odd  values 
of  N  approach  ^  from  below,  while  the  sums  for  even  values  of  N  approach  it  from 
above.  For  negative  integer  values  of  k ,  the  partial  sums  approach  —  \  instead. 
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for  even  and  odd  values  of  N  up  to  a  million.  Assuming  S(N)  — >  as  N  — >  oo  for  all  k  >  0,  then 


ui  =  V2Glu o  «  1.353wo 


(37) 


where  the  Catalan  constant  is  G  =  (3{ 2)  =  1  —  1/9  +  1/25  —  ...  ~  0.91597.  This  approximate 
frequency  is  8%  larger  than  the  exact  value,  which  is  an  improvement  over  the  standard  harmonic 
balance  estimate  obtained  as  follows.  Equation  (1)  can  be  rewritten  for  the  inverse-distance  force 
law  as 


dt2 


=  -A2 


(38) 


Substitute  x  =  A  cos  ujt  to  get 

|w2(l  +  cos2wf)  =  u>q,  (39) 

using  the  double-angle  formula  for  cosine,  and  balance  the  dc  terms  to  find 


u>  =  V2u)o  «  1.414wo 


(40) 


which  is  13%  larger  than  the  exact  value  in  Eq.  (31). 


4.4  The  large-amplitude  simple  pendulum 

A  point  mass  m  is  suspended  from  an  ideal  string  of  length  L  whose  angular  displacement  is  9 
relative  to  the  vertical.  Defining  Wg  =  g/L ,  the  restoring  torque  divided  by  the  moment  of  inertia 
of  the  bob  (both  about  the  pivot)  is  the  trigonometric  function  f{6)  =  WgSin6f  Inserting  this 
torque  into  Eq.  (19)  with  x  replaced  by  9 ,  the  exact  angular  frequency  is  found  to  be 


k^exact 


2-i/2 


TTUIq 


l-A  dO 

JO  y/cOsW^COsA 


(41) 


which  can  be  expressed  in  terms  of  an  elliptic  integral.  Evaluating  it  for  an  angular  amplitude  of 
A  =  7t/2  gives 


k^exact 


2-7T3/2 

r2(i/4)w° 


0.8472cjq 


(42) 


where  T  is  the  gamma  function.  (The  reciprocal  of  this  frequency  implies  that  the  period  of 
a  pendulum  swinging  between  90°  and  +90°  is  18%  longer  than  it  would  be  at  small  angular 
amplitudes.) 

Equation  (3)  becomes 


27m  (— l)n+1  sin  A  2 
n27r2  -  A2  W°' 


(43) 


Equation  (13)  then  predicts 


to 


2 


4wg  sin  A 
nA 


OO 


E 


{2k  +  l)2 


E 


(— l)n+1n  J2fc+i(mr) 
n2  —  A2/tt2 


(44) 


Expression  19  on  page  679  of  Prudnikov,  Brychkov,  and  Marichev  [8]  is 


E 


(— l)n+1n  J2k+\{nx) 
n2  —  a 2 


7 r 
2 


csc(a7r)  J2fe+i(aa;). 


(45) 
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Evaluating  it  at  a  =  A/- tt  and  x  =  tt,  and  substituting  the  result  into  Eq.  (44),  one  obtains 


2  _  (—  l)fc'/2fc+l(^) 

"  =  A  h  (2fc+1)2 


(46) 


which  reproduces  Eq.  (25)  of  Ref.  [9].  Retaining  only  the  first  term  in  this  series,  the  period  of 
oscillation  is  found  to  be 


TT77H  (47) 


where  To  =  2tt^J L/ g  is  the  familiar  period  for  small  amplitudes.  This  result  agrees  with  the 
standard  harmonic  balance  method,  Eq.  (34)  of  Ref.  [10].  Equation  (46)  can  be  thought  of  as  its 
extension  to  all  overtones.  For  the  large-amplitude  case  A  =  7r/2,  it  becomes 


(48) 


There  is  no  closed  form  for  the  sum  of  this  Neumann  series.  However,  using  the  standard  series 
expansion  for  a  Bessel  function  of  positive  integer  order,  one  immediately  sees  that  ^fc+iX71"/^)  is 
positive  and  monotonically  decreasing  in  value  as  k  increases.  Thus  Eq.  (48)  is  an  alternating  series 
of  terms  of  decreasing  magnitude.  It  follows  that  one  can  cut  off  the  series  at  an  upper  limit  of 
k  =  N  and  get  an  approximate  sum  with  an  error  of  less  than  the  first  omitted  term.  Choosing 
N  =  2,  the  error  is  already  less  than  J7(7t/2)/49  w  7x  10-7,  in  which  case 


(49) 


w  «  0.8438u;o 


in  good  agreement  with  Eq.  (42). 

4.5  The  Duffing  oscillator  with  zero  linear  term 

For  an  anharmonic  oscillator  having  restoring  force  fix)  =  ax3,  define  ui o  =  A\Ja.  Using  the 
change  of  variables  x  =  Au,  Eq.  (19)  gives  the  same  frequency  as  Eq.  (42). 

On  the  other  hand,  Eq.  (3)  becomes 


(50) 


so  that  Eq.  (13)  predicts 


(51) 


Using  Eqs.  (55),  (56),  (58),  (59),  and  (65)  in  Appendix  A,  this  prediction  becomes 


(52) 


which  is  4%  larger  than  the  exact  value  in  Eq.  (42).  Interestingly,  the  predicted  frequency  becomes 
-\/3  wo/2  if  one  retains  only  the  J\  terms  (corresponding  to  k  =  0)  in  the  sums  in  Eq.  (51),  reducing 
the  error  to  2%. 
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5  Conclusions 

In  this  paper,  a  double  Fourier  harmonic  balance  method  has  been  developed,  based  on  the 
expansion  of  the  force  in  a  Fourier  sine  series.  By  applying  this  new  method  to  a  number  of 
previously  studied  nonlinear  oscillators,  approximate  frequencies  have  been  obtained  that  are  often 
better  than  those  found  using  other  techniques  [11].  The  method  works  whenever  the  Fourier 
coefficients  in  Eq.  (3)  and  the  sums  of  the  Bessel  series  in  Eq.  (13)  can  be  calculated,  although 
a  good  approximation  can  often  be  obtained  by  retaining  only  the  J\  term  in  the  latter  equation. 
As  additional  examples  beyond  the  cases  treated  here,  the  sublinear  oscillator  f(x)  oc  x 1/3  and  the 
general  power-law  oscillator  ±xn  with  n  >  3  can  be  expressed  in  terms  of  the  incomplete  gamma 
function.  The  method  should  also  be  applicable  to  damped  oscillators,  by  including  the  full  (not  just 
the  sine  and  odd  cosine)  expansions  in  Eqs.  (2)  and  (5).  As  a  consequence,  this  new  version  of  the 
harmonic  balance  method  is  pedagogically  useful  in  analyzing  the  behavior  of  nonlinear  oscillating 
systems.  The  derivations  are  accessible  to  undergraduate  students,  in  contrast  to  say  the  homotopy 
method  [12,  13],  especially  because  of  the  ready  availability  of  Fourier  analysis  and  Bessel  functions 
in  computer  algebraic  software. 

Appendix  A 

In  this  appendix,  the  sum 


(53) 


is  evaluated,  where  k  is  any  positive  odd  integer,  and  p  equals  either  1  or  3.  There  are  three  cases 
of  interest  as  follows. 

Case  1:  k  =  p 

Equation  (2)  on  page  635  of  Watson  [14]  can  be  rearranged  as 


‘it—  x 


(54) 


where  m  is  any  positive  integer  and  0  <  x  <  tt.  In  particular,  the  two  needed  values  are 


(55) 


n= 1 


and 


n= 1 


(56) 


Case  2:  k  > 


Equation  (3)  on  page  636  of  Watson  [1  ]  can  be  rearranged  as 


(57) 
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where  m  >  2s  with  m  and  s  both  positive  integers,  and  0  <  x  <  n.  In  particular,  one  set  of  values 
of  interest  is 

S(k,  1)  =  ^1)n+^k<yn^  =  0  for  *  =  3, 5, 7, ...  (58) 

n= 1 

as  obtained  by  substituting  {to  =  3  &  s  =  1},  {to  =  5  &  s  =  2},  {to  =  7  &  s  =  3},  . . .  into  Eq.  (57). 
The  other  needed  set  of  values  is 

S(k,  3)  =  £  =  0  for  k  =  5, 7, 9, ...  (59) 

n—1  '  7 

as  obtained  by  substituting  {to  =  5  &  s  =  1},  {to  =  7  &  s  =  2},  {to  =  9  &  s  =  3},  . . .  into  Eq.  (57). 
Case  3:  k  =  1  &  p  =  3 


The  standard  Bessel  recursion  relation  can  be  written  as 

z  z 

Js (z)  =  —Js-l(z)  +  —  Js+1(z). 


(60) 


Multiply  this  equation  through  by  (— l)n+1z  m,  let  z  =  nx,  and  then  sum  over  n  to  get 


y,  (-1  )n+lJ8(nx)  =  l_  y.  (-l)n+1Js-i(nx)  1  y>  (-l)n+1  Js+1  (nx) 

( TlX )m  9.q  (n/r\rn~  1  9<?  Z-^d 


(nx) 


(nx) 


m—  1 


n—1  v  7  n—1  v  7  n—1 

Next  substitute  s  =  1,  to  =  3,  and  x  =  7r  so  that  this  result  specifically  becomes 

^  (— 1)"+1J1(n7r)  =  1^  (— 1)”+1  Jo(n7r)  1^,  _(-l)"+1  J2(n7r) 
(mr)3  2  (mr)2  2 


(tut)2 


(61) 


(62) 


According  to  Eq.  (54),  the  third  sum  is  5(2,2)  =  1/16.  To  get  the  other  sum,  use  the  Schlomilch 
series  [15] 


(— 1  )n+1  J0(nx)  7T2  x2 

r)2 


n—1 


12 


Divide  through  by  x2  and  then  let  x  =  ir  to  obtain 


5(0,2)  =  ^ 


n—1 


(-l)ra+1  J0(?r7r) 
(mr)2 


1  1 

12  ~~  8 


1 

'24' 


(63) 


(64) 


Finally  substitute  these  numerical  values  of  5(2,  2)  and  5(0,  2)  into  the  right-hand  side  of  Eq.  (62) 
to  get 


5(1.3) -f:  (-i)ri$(”)q 

'  (mrr  2 

n—1  v  7 


-l  l 
24  +  16 


1 

96 


(65) 


which  happens  to  equal  5(3,3)  from  Eq.  (56). 
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Appendix  B 


The  Schlomilch  series  for  x  is  [15] 


=  —  —  2  V  (66) 

A  n 2  v  ' 


odd  n—  1 


Take  its  derivative  with  respect  to  x  to  get 


l  =  -2  E 


However 


and  thus 


1  cLJo(nx) 
n  d(nx ) 

odd  n—  1  '  ' 

dJp  jz) 
dz 


=  ~Ji  0) 


(67) 


(68) 


E  ^  =  2-  (69) 

odd  n=  1 

Remarkably,  the  sum  is  independent  of  x  over  the  range  from  0  to  n.  In  particular,  it  can  be 
evaluated  at  x  =  n  to  get  the  result  needed  in  Eq.  (26). 


Appendix  C 


Expression  (i)  of  Petkovic  [16]  with  k  — >  k  +  1  and  m  =  1  becomes 

J2k+i(nx)  1 


E 

n—l 


n  2k  +  1 


for  integers  k  >  1  and  0  <  x  <  2tt.  Evaluate  Eq.  (70)  at  x  =  tt  to  obtain 

J2fc+i(wr)  1 


E 


n  2k  +  1 


Next  put  n  =  p/2  into  Eq.  (70)  to  get 

OO 

E 


J2k+i{px/2)  1 


even  p—2 


and  substitute  x  =  2-k  to  find 


E 

even  n— 2 


p  2(2k  +  1)  ’ 


J2k+i(nir)  1 


n  2(2k  +  1) 


Finally,  subtract  Eq.  (73)  from  (71)  to  establish  that 

J2fc+i(n7r)  1 


E 

odd  n—l 


n  2(2k  +  1) ' 


(70) 

(71) 

(72) 

(73) 

(74) 


Remarkably,  Eq.  (74)  is  valid  for  k  =  0,  as  Eq.  (69)  shows,  even  though  neither  Eq.  (71)  nor  (73) 
is  valid  for  k  =  0  (as  one  can  verify  by  numerically  approximating  the  latter  two  sums  for  k  =  0). 
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